<html>
<link href="../enzo.css" rel="stylesheet" type="text/css">
  <head>
    <title>Enzo Algorithms</title>
  </head>
<body> 
    <h1>Enzo Algorithms</h1>

    <p>This section provides a very short overview of the algorithms used by the Enzo code.
      References to texts and journal articles providing a more complete discussion of the
      algorithms are included at the end of this page for the interested reader.  As of this
      writing (January 2004) a formal Enzo method paper has not been published, but is in 
      preparation.  Much of the text and images on this page have been taken from one of the 
      <a href="http://cosmos.ucsd.edu/">Laboratory for Computational Astrophysics</a> contributions 
      to the <a href="http://flash.uchicago.edu/amr2003/">2003 AMR conference</a>.<b>[1]</b>      
    </p>

<p>
Enzo is written in a mixture of C++ and Fortran 77.  High-level functions and data structures
are implemented in C++ and computationally intensive lower-level functions are written in Fortran.
Enzo is parallelized using the <a href="http://www-unix.mcs.anl.gov/mpi/">MPI</a> message-passing
library and uses the <a href="http://hdf.ncsa.uiuc.edu/HDF5/">HDF 5</a> data format to write out
data and restart files in a platform-independent format.
</p>

<h3>Adaptive Mesh Refinement</h3>
<p>
Enzo allows hydrodynamics in 1, 2 and 3 dimensions using the structured adaptive mesh refinement
(SAMR) technique developed by Berger and Collela.  <b>[2]</b>  The code allows arbitrary integer ratios
of parent and child grid resolution and mesh refinement based on a variety of criteria, including baryon
and dark matter overdensity or slope, the existence of shocks, Jeans length, and cell cooling time.  The
code can also have fixed static nested subgrids, allowing higher initial resolution in a subvolume of 
the simulation.  Refinement can occur anywhere within the simulation volume or in a user-specified 
subvolume.</p>

<p>
The AMR grid patches are the primary data structure in \enzo.  Each individual 
patch is treated as an individual object, and can contain both field variables
and particle data.  Individual patches are organized into a dynamic distributed
AMR mesh hierarchy using arrays of linked lists to pointers to grid objects.
The code uses a simple dynamic load-balancing scheme to distribute the workload 
within each level of the AMR hierarchy evenly across all processors.
</p>
<p>
Although each processor stores the entire distributed AMR hierarchy, not all
processors contain all grid data.  A grid is a <i>real grid</i> on a 
particular processor if its data is allocated to that processor, and a 
<i>ghost grid</i> if its data is allocated on a different processor.  Each
grid is a real grid on exactly one processor, and a ghost grid on all others.
When communication is necessary, MPI is used to transfer the mesh or particle
data between processors.  The tree structure of a small illustrative 2D AMR
hierachy - six total grids in a three level hierarchy distributed across
two processors - is shown on the left in Figure 1.
</p>
<p>
<center><img src="amr_hierarchy.jpg"><br>
<b>Figure 1.  </b>  Real and ghost grids in a hierarchy; real and ghost zones in a grid.</center>
</p>

<p>
Each data field on a real grid is an array of zones with dimensionality
equal to that of the simulation (typically 3D in cosmological structure
formation).  Zones are partitioned into a core block of <i>real zones</i>
and a surrounding layer of <i>ghost zones</i>.  Real zones are used to
store the data field values, and ghost zones are used to temporarily store
values from surrounding areas, ie, neighboring grids, parent grids or external
boundary conditions, when required for updating real zones.
The ghost zone layer is three zones deep in order to accomodate the computational
stencil in the hydrodynamics solver (See below), as indicated
in the right panel in Figure 1.  These ghost
zones can lead to significant computational and storage overhead, especially
for the smaller grid patches that are typically found in the deeper levels of an 
AMR grid hierarchy.
</p>
<p>
For more information on Enzo implementation and data structures, see 
references <b>[3], [4], [5]</b> and <b>[6]</b>
</p>

<h3>Dark Matter Dynamics</h3>
<p>

The dynamics of large-scale structures are dominated by dark matter,
which accounts for approximately 85% of the matter in the universe but can only
influence baryons via gravitational interaction. 
There are many other astrophysical situations where gravitational 
physics is important as well, such as galaxy collisions, where the stars in
the two galaxies tend to interact in a collisionless way.
</p>
<p>
Enzo uses the Particle-Mesh N-body method to calculate collisionless particle dynamics.
This method follows trajectories of a representative 
sample of individual particles and is much more efficient than a direct 
solution of the Boltzmann equation in most astrophysical situations.
The particle trajectories are controlled by a simple set of coupled equations
(for simplicity, we omit cosmological terms):
<p>
<img src="graveqtn1.jpg">
</p>
and
<p>
<img src="graveqtn2.jpg">
</p>
<p>
Where x<sub>p</sub> and v<sub>p</sub> are
 the particle position and
velocity vectors, respectively, and the term on the right-hand side of 
the second equation is the gravitational force term.  The solution to this
can be found by solving the elliptic Poisson's equation:</p>

<p>
<img src="graveqtn3.jpg">
</p>
<p>
where &rho; is the density of both the collisional fluid (baryon gas)
and the collisionless fluid (particles).
</p>
<p>
These equations are finite-differenced and for simplicity are 
solved with the same timestep as the equations of hydrodynamics.
The dark matter particles are sampled onto the grids using the 
triangular-shaped cloud (TSC) interpolation technique to form a
spatially discretized density field (analogous to the baryon
densities used to calculate the equations of hydrodynamics)
and the elliptical equation is solved using FFTs on the
triply periodic root grid and multigrid relaxation on the subgrids.
Once the forces have been computed on the mesh, they are interpolated
to the particle positions where they are used to update their 
velocities.
</p>


<h3>Hydrodynamics</h3>
<p>
The primary hydrodynamic method used in Enzo is based on the 
piecewise parabolic method (PPM) of Woodward & Colella <b>[7]</b>
which has been significantly modified for the study of cosmology.  The
modifications and several tests are described in much more detail in <b>[8]</b>,
and we recommend that the interested reader look there.
</p>
<p>
PPM is a higher-order-accurate version of Godunov's method with 
third-order-accurate piecewise parabolic monotolic interpolation 
and a nonlinear Riemann solver for shock capturing.  It does an excellent
job capturing strong shocks and outflows.  Multidimensional schemes are 
built up by directional splitting, and produce a method that is formally 
second-order-accurate in space and time and explicitly conserves energy, 
momentum and mass flux.  The conservation laws for fluid mass, momentum
and energy density are written in comoving coordinates for a
Friedman-Robertson-Walker spacetime.  Both the conservation laws and
Riemann solver are modified to include gravity, which is calculated as
discussed above.
</p>
<p>

There are many situations in astrophysics, such as the bulk hypersonic
motion of gas, where the kinetic energy of a fluid can dominate its internal 
energy by many orders of magnitude.  In these situations, limitations on
machine precision can cause significant inaccuracy in the calculation of 
pressures and temperatures in the baryon gas.  In order to address this 
issues, Enzo solves both the internal gas energy equation and the total 
energy equation everywhere on each grid, at all times.  This 
<i>dual energy formalism</i> ensures that the method yields the correct
entropy jump at strong shocks and also yields accurate pressures and 
temperatures in cosmological hypersonic flows.  See reference <b>[8]</b> for
more information about the dual energy formalism.
</p>
<p>
As a check on our primary hydrodynamic method, we also include an
implementation of the hydro algorithm used in the Zeus astrophysical
code.  <b>[9], [10]</b>  This staggered grid, finite difference
method uses artificial viscosity
as a shock-capturing technique and is formally first-order-accurate when
using variable timesteps (as is common in structure formation simulations), 
and is not the preferred method in the Enzo code.
</p>

<h3>Cooling/Heating</h3>
<p>
The cooling and heating of gas is extremely
important in astrophysical situations.  To this extent, 
two radiative cooling models and several uniform ultraviolet background 
models have been implemented in an easily extensible framework.
</p>
<p>
The simpler of the two radiative cooling models assumes that all
species in the baryonic gas are in equilibrium and calculates cooling 
rates directly from a cooling curve assuming Z = 0.3 Z<sub>o</sub>.  The 
second routine, developed by Abel, Zhang, Anninos & Norman <b>[11]</b>,
assumes that the gas has primordial abundances (ie, a gas which is composed of
hydrogen and helium, and unpolluted by metals), and solves
a reaction network of 28 equations which includes collisional and radiative processes
for 9 seperate species (H, H<sup>+</sup>, He, He<sup>+</sup>, He<sup>++</sup>, H<sup>-</sup>,
 H<sub>2</sub><sup>+</sup>, H<sub>2</sub>, and e<sup>-</sup>).  In 
order to increase the speed of the calculation, this method takes the reactions with 
the shortest time scales (those involving H<sup>-</sup> and H<sub>2</sub><sup>+</sup>) 
and decouples them from
the rest of the reaction network and imposes equilibrium concentrations, which is 
highly accurate for cosmological processes.  See <b>[11]</b> and 
<b>[12]</b> for more information.

The vast majority of the volume of the present-day universe is occupied by
low-density gas which has been ionized by ultraviolet radiation from quasars, 
stars and other sources.  This low density gas, collectively referred to as the
``Lyman-&alpha; Forest'' because it is primarily observed as a dense collection of
absorption lines in spectra from distant quasars (highly luminous extragalactic
objects), is useful because it can be used to determine several cosmological
parameters and also as a tool for studying the formation and evolution of 
structure in the universe (see <b>[13]</b> for more information).  The 
spectrum of the ultraviolet radiation background plays an important part
in determining the ionization properties of the Lyman-&alpha; forest, so 
it is very important to model this correctly.  To this end, we have
implemented several models for uniform ultraviolet background radiation
based upon the models of Haardt & Madau <b>[14]</b>.
</p>
<h3>Star Formation/Feedback</h3>
<p>
One of the most important processes when studying the formation and
evolution of galaxies (and to a lesser extent, groups and clusters of 
galaxies and the gas surrounding them) is the formation and feedback
of stars.  We use a heuristic prescription similar to that of 
Cen & Ostriker <b>[15]</b> to convert gas which is rapidly
cooling and increasing in density into star ``particles'' which
represent an ensemble of stars.  These particles then evolve 
collisionlessly while returning metals and thermal energy back
into the gas in which they formed via hot, metal-enriched winds.  

</p>
<h3>Parallelization in Enzo</h3>
<p>
Enzo uses a grid-based parallelization scheme for load balancing.  The root grid is 
partitioned up into N pieces (where N is the number of processors), and each processor
is given a piece of the root grid, which it keeps for the duration of the simulation run.
Subgrids are treated as independent objects and are distributed to the processors such that
each level of grids is load-balanced across all processors.  Boundary fluxes between neighboring
grid patches and parent and children grids are passed back and forth using MPI commands.</p>
<p>
The one portion of the code that is parallelized differently is the root grid gravity solver.  As
discussed above, the gravitational potential on the root grid is solved using a fourier
transform method, which requires its own message-passing routines.  The three-dimensional total
density field (composed of the dark matter plus baryon density on the root grid) is decomposed 
into two-dimensional slabs (requiring one set of messages), which are then fourier transformed.  
The slabs are then transposed along another axis (requiring a second set of messages to be passed)
and transformed again, and a third set of messages is required in order to obtain the original
block decomposition.  This is unavoidable when using a fourier transform scheme, and as a result
the speed of the root grid gravity solver is very sensitive to the speed of the communication
network on the platform that Enzo is being run on.
</p>

<h3>Initial Conditions Generator</h3>
<p>
A somewhat detailed description of the method Enzo uses to create initial conditions can be 
downloaded as a <a href="makeics.ps">postscript</a> or <a href="makeics.pdf">PDF</a> document.
To summarize:  Dark matter particles and baryon densities are laid out on a uniform Cartesian 
grid.  Given a user-specified power spectrum P(k), the linear density fluctuation field is calculated at
some initial time (typically z = 100 for high-resolution/small box simulations) by using P(k) to obtain
the density fluctuations in k-space on a uniform Cartesian grid.  P(k) is sampled discretely at each
grid point, with the density fluctuations having a random complex phase and amplitude.  The amplitude
is generated such that the distribution of amplitudes is Gaussian.  This cube is then fourier transformed
to give physical density fluctuations.  Particle positions and velocities and baryon velocities are
calculated using the Zel'Dovich approximate.  See the document above, or read Bertschinger 1998 <b>[16]</b>
for more information.
</p>


<h3>References</h3>
<p>
<b>[1]</b> B. W. O'Shea et al.  "Introducing Enzo, an AMR Cosmology
 Application." To be published in <u>Adaptive Mesh Refinement - Theory And Applications</u>,
 the proceedings from the 2003 University of Chicago AMR Workshop<br>
<b>[2]</b> M. J. Berger and P. Colella.  "Local adaptive mesh refinement for shock hydrodynamics," 
<i>J. Comp. Phys</i>, 82:64-84, 1989<br>
<b>[3]</b> G. L. Bryan.  "Fluids in the universe: Adaptive mesh in Cosmology."  
<i>Computing in Science and Engineering</i>, 1:2, 1999<br>
<b>[4]</b> G. L. Bryan and M. L. Norman.  "A hybrid AMR application for cosmology and astrophysics." 
In <i>Workshop on Structured Adaptive Mesh Refinement Grid Methods"</i>, p. 165.  
IMA Volumes in Mathematics #117, 2000<br>
<b>[5]</b> G. L. Bryan and M. L. Norman.  In D.A. Clarke and M. Fall, editors, 
<u>Computational Astrophyiscs: 12th Kingston Meeting on Theoretical Astrophysics, 
proceedings of a meeting held in Halifax; Nova Scotia; Canada Oct. 17-19, 1996.</u> 
ASP Conference Series #123, 1997<br>
<b>[6]</b> M. L. Norman and G. L. Bryan. "Cosmological Adaptive Mesh Refinement." 
In Kohji Tomisaka, Shoken M. Miyama and Tomoyuki Hanawa, editors, 
<u>Numerical Astrophysics: Proceedings of the International Conference on Numerical Astrophysics 1998</u>, 
p. 19.  Kluwer Academics, 1999<br>
<b>[7]</b> P. R. Woodward and P. Colella.  "A piecewise parabolic method for gas dynamical simulations," 
<i>J. Comp. Phys</i>, 54:174, 1984<br>
<b>[8]</b> G. L. Bryan, M. L. Norman, J. M. Stone, R. Cen and J. P. Ostriker. "A piecewise 
parabolic method for cosmological hydrodynamics," 
<i>Comp. Phys. Comm.</i>, 89:149, 1995<br>
<b>[9]</b> J. M. Stone and M. L. Norman.  "Zeus-2D: A radiation magnetohydrodynamics code for 
astrophysical flows in two space dimensions.  I.  The hydrodynamics algorithms and tests." 
<i>The Astrophysical Journal</i>, 80:753, 1992<br>
<b>[10]</b> J. M. Stone and M. L. Norman.  "Zeus-2D: A radiation magnetohydrodynamics code for 
astrophysical flows in two space dimensions.  II. The magnetohydrodynamic algorithms and tests."
 <i>The Astrophysical Journal</i>, 80:791, 1992 <br>
<b>[11]</b> T. Abel, P. Anninos, Y. Zhang and M.L. Norman.  "Modeling primordial gas in numerical cosmology." 
<i>New Astronomy</i>, 2:181-207, August 1997<br>
<b>[12]</b> P. Anninos, Y. Zhang, T. Abel and M.L. Norman.  "Cosmological hydrodynamics with multispecies chemistry 
and nonequilibrium ionization and cooling."  <i>New Astronomy</i>, 2:209-224, August 1997<br>
<b>[13]</b> M. Rauch.  "The Lyman Alpha Forest in the Spectra of QSOs." <i>Annual Review of Astronomy and 
Astrophysics</i>, 36:267-316, 1998<br>
<b>[14]</b> F. Haardt and P. Madau.  "Radiative Transfer in a Clumpy Universe, II. The Ultraviolet Extragalactic 
Background."  <i>The Astrophysical Journal</i>, 461:20, 1996<br>
<b>[15]</b> R. Cen and J.P. Ostriker. "Galaxy formation and physical bias." <i>The Astrophysical Journal Letters</i>, 399:L13, 1992<br>
<b>[16]</b> E. Bertschinger.  "Computer Simulations in Cosmology." <i>Annual Review of Astronomy and Astrophysics</i>, 36:599<br>
</p>

<p>&nbsp;</p>
<p>
<a href="index.html">Previous - Index</a><br>
<a href="ics_top.html">Next - Generating Initial conditions</a><br>
</p>


<p>&nbsp;</p>
<p>
<a href="../index.html">Go to the Enzo home page</a>
</p>
<hr WIDTH="100%">
<center>&copy; 2004 &nbsp; <a href="http://cosmos.ucsd.edu">Laboratory for Computational Astrophysics</a><br></center>
<center>last modified February 2004<br>
by <a href="mailto:bwoshea (AT) lanl.gov">B.W. O'Shea</a></center>


</body>
</html>
